subroutine getpsivel
use mod_2dflu
implicit none
integer :: i1,i2,k1,k2,ireal,iimag,ksqr
double precision :: rk2,kx,ky
double complex :: oo,uu_x,uu_y,psi_complex
!! get psi and velocity from the omega
factor = 2.0d0*pi/length

do i2 = 1,n2
  k2 = (i2-1) - n2*(i2/(n1hf+1))	
  do i1 =1,n1hf
    k1 = i1 -1
    ksqr = k1*k1 + k2*k2 
    rk2 = factor*factor*dfloat(ksqr)
    ireal = 2*i1 -1
    iimag = 2*i1
    kx = factor*dfloat(k1)
    ky = factor*dfloat(k2)
    if((ksqr.gt.kasqr).or.(ksqr.eq.0))then
      omega(ireal,i2) = 0.0d0
      omega(iimag,i2) = 0.0d0
      psi(ireal,i2) = 0.0d0
      psi(iimag,i2) = 0.0d0
      ukx(ireal,i2) = 0.0d0
      ukx(iimag,i2) = 0.0d0
      uky(ireal,i2) = 0.0d0
      uky(iimag,i2) = 0.0d0
    else
      oo=dcmplx(omega(ireal,i2),omega(iimag,i2))
!
! \lap \psi = -\omega, -k^2 \hat{\psi} = -\hat{\omega},\hat{\psi} = \hat{\omega}/k^2
!
      psi_complex=oo/rk2
!
! u_x = \partial_y \psi, \hat{u}_x = i ky \hat{psi}
!
      uu_x = zi*ky*psi_complex
!
! u_y = -\partial_x \psi, \hat{u}_x = -i kx \hat{psi}
!
      uu_y = -zi*kx*psi_complex
!
      psi(ireal,i2) = real(psi_complex);psi(iimag,i2) = aimag(psi_complex)
      ukx(ireal,i2) = real(uu_x);       ukx(iimag,i2) = aimag(uu_x) 
      uky(ireal,i2) = real(uu_y);       uky(iimag,i2) = aimag(uu_y) 
    endif
  enddo
enddo
!! --------------------- 
end subroutine getpsivel
